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Diamond is studied by path integral molecular dynamics simulations of the atomic nuclei in com- 
bination with a tight-binding Hamiltonian to describe its electronic structure and total energy. This 
approach allows us to quantify the influence of quantum zero-point vibrations and finite tempera- 
tures on both the electronic and vibrational properties of diamond. The electron-phonon coupling 
mediated by the zero-point vibration reduces the direct electronic gap of diamond by 10%. The 
calculated decrease of the direct gap with temperature shows good agreement with the experimental 
data available up to 700 K. Anharmonic vibrational frequencies of the crystal have been obtained 
from a linear-response approach based on the path integral formalism. In particular, the temper- 
ature dependence of the zone-center optical phonon has been derived from the simulations. The 
anharmonicity of the interatomic potential produces a red shift of this phonon frequency. At tem- 
peratures above 500 K, this shift is overestimated in comparison to available experimental data. 
The predicted temperature shift of the elastic constant C44 displays reasonable agreement with the 
available experimental results. 

PACS numbers: 63.20.Kr, 81.05.Uw, 63.20.Ry, 05.30.-d 
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I. INTRODUCTION 

Tctrahedral semiconductors such as silicon, germa- 
nium or diamond, have served as model materials to 
study the electronic and vibrational properties of crys- 
tals. In particular, the effects of the lattice vibrations 
on the electronic properties, through the mechanism of 
electron-phonon coupling, have been experimentally in- 
vestigated by measuring the temperature dependence of 
their optical excitation spectrai Furthermore, the ob- 
servation of the dependence of such spectra with isotopic 
mass has provided detailed information on the electronic 
properties of these solidsi^ Besides the electron-phonon 
interaction, the anharmonicity of lattice vibrations has 
been observed by the dependence of the phonon fre- 
quencies and line-widths with temperature and isotopic 
composition.^ 

From a theoretical point of view, in spite of the impres- 
sive progress of ab initio methods for the investigation of 
the electronic structure of solids, the atomic nuclei are 
usually considered either as fixed in their crystallographic 
positions or by approximating their dynamics with clas- 
sical mechanics.'^ Thus, the effects of the electron-phonon 
interactions on electronic properties and the effect of 
zero-point anharmonicity in the vibrational properties of 
the lattice are usually neglected in these calculations. As 
it has been pointed out in Ref. Q, these effects may be 
even larger that the error assumed in the electronic ab 
initio calculations. 

The electron-phonon interaction in tetrahedral semi- 
conductors has been studied theoretically by perturba- 
tion theory.A^ Within this approach, the reduction of the 



direct electronic gap due to zero-point vibrations of the 
lattice phonons is predicted to be of 0.62 eV in diamond 
and of 0.06 eV in germanium. Both energy shifts repre- 
sent roughly a 10% fraction of the corresponding energy 
gaps, i.e., they are so large that a quantitative descrip- 
tion of the electronic structure can not be expected by 
a theory that neglects such effects. Anharmonic shifts 
of the phonon modes of diamond and silicon have been 
determined by combining density-functional perturba- 
tion theory with a frozen-phonon approach,^ and the re- 
sults show good agreement with the experimental data 
available from first-order Raman spectrai^ A review of 
the current status of lattice-dynamical calculations using 
density-functional perturbation theory can be found in 
Ref. 13 

The path integral (PI) formulation of statistical me- 
chanics offers an alternative way to study finite tempera- 
ture properties that are related to the quantum nature of 
the atomic nucleiii2iii Thus, the combination of the path 
integral formulation with electronic structure methods is 
an interesting alternative to perturbational approaches 
for the study of electronic and vibrational properties of 
solids. An advantage of this approach is that both the 
electrons and the atomic nuclei are treated quantum me- 
chanically in the framework of the Born-Oppenheimer 
(BO) approximation, so that anharmonic and tempera- 
ture effects can be evaluated for both vibrational proper- 
ties and the electronic structure. This unified scheme has 
been applied so far to the study of solids and molecules 
containing light atomsi2J^JiiiS*iSiiLi2iia 

In this paper we present a path integral molecular dy- 
namics study of diamond at temperatures between 100 
and 1200 K. The electronic structure has been treated 



by a non-orthogonal tight-binding (TB) Hamiltonian as 
a reasonable compromise to reduce the computational 
cost of deriving the BO surface for the nuclear dynam- 
ics. In particular, we are interested in the investigation of 
electronic properties that are determined by the electron- 
phonon coupling, as the dependence of the electronic gap 
with temperature and isotopic mass. Also vibrational 
properties that depend on the anharmonicity of the in- 
teratomic potential will be studied with the help of a 
linear response approach recently developed within the 
path integral formalism,^'' PI simulations of diamond us- 
ing effective interaction potentials have been carried out 
earlier to study structural and thermodynamic properties 
of this material.-- 

This paper is organized as follows. In Sec. II, we de- 
scribe the computational method and the models em- 
ployed in our simulations. Our results are presented and 
discussed in Sec. Ill, dealing with the direct electronic 
gap, the vibrational energy of the solid, and the temper- 
ature dependence of the frequency of the optical phonon 
at the center of the Brillouin zone (BZ) and the elastic 
constant C44 of diamond. In Sec. IV, we present the main 
conclusions of the paper. 



II. COMPUTATIONAL METHOD 
A. Simulation details 

The formalism employed here for the quantum treat- 
ment of electrons and nuclei is based on the combina- 
tion of the path integral formulation, to derive proper- 
ties of the atomic nuclei in thermal equilibrium, with an 
electronic Hamiltonian to derive the BO energy surface, 
E 30(^)1 a-s a function of the nuclear configuration R. 
The path integral and electronic structure parts of the 
resulting algorithm appear as independent blocks, since 
the only electronic result required for the path integral 
simulation is the value of the function _Bbo(R), and pos- 
sibly its derivatives with respect to ionic positions. Thus, 
the combination of path integrals with any chosen elec- 
tronic Hamiltonian is straightforward. For the present in- 
vestigation of diamond we have chosen an efhcient tight- 
binding one-electron effective Hamiltonian, based on den- 
sity functional (DF) calculations>22i The use of a simpli- 
fied electronic Hamiltonian is a reasonable compromise to 
explore the efficiency and capability of this unified for- 
malism for the evaluation of electronic and vibrational 
properties of solids at finite temperatures. The imple- 
mentation of density functional or Hartree-Fock based 
self-consistent methods is left for future development. 
The capability of tight-binding methods to simulate dif- 
ferent properties of solids and molecules has been re- 
viewed by Goringe et alM' 

The computational advantage of using the path- 
integral formulation of statistical mechanics is formulated 
by the so-called " quantum-classical" isomorphism. Thus, 
this method exploits the fact that the partition function 
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FIG. 1: Vibrational energy of diamond as a function of the 
time step. At, employed in the PI MD algorithm. The results 
are derived at 300 K with a Trotter number L — 20. The 
continuous line is a quadratic fit to the simulation results. 



of a quantum system is formally equivalent to that of a 
classical one, obtained by replacing each quantum parti- 
cle (here, atomic nucleus) by a ring polymer consisting 
of L "beads", connected by harmonic springs >i2iii*2iiS^ 
In many-body problems, the configuration space of the 
classical isomorph is usually sampled by Monte Carlo or 
molecular dynamics (MD) techniques. Here, we have em- 
ployed the PI MD method, which has been found to re- 
quire less computer time resources in our problem. Ef- 
fective algorithms to perform PI MD simulations in the 
canonical NVT ensemble have been described in detail 
by Martyna et al?^ and TuckermaniSi All calculations 
presented here were carried out in the canonical ensem- 
ble, using originally developed MD software, which en- 
ables efficient PI MD simulations on parallel supercom- 
puters. 

Simulations were carried out on a 2 x 2 x 2 super- 
cell of the diamond face-centered cubic cell with periodic 
boundary conditions, containing iV = 64 C atoms. The 
atomic mass of carbon was set to 12 amu. The conver- 
gence of the internal energy has been checked for some 
selected atomic configurations, by considering sets of 1, 
4 and 32 k points in the BZ of the simulation supercell. 
The main effect of the k point sampling is found to be 
a constant shift of the internal energy. This systematic 
error is largely reduced in the calculation of properties 
obtained as energy differences (e.g., energy shifts as a 
function of temperature). For this reason, we have cho- 
sen to use only the F point for the sampling of the BZ 
of the simulation supercell. A set of 4 k points increases 
the computer time by a factor of 10 with respect to the F 
point sampling, without significant changes of the results 
presented here. 

The simulation-cell parameter employed in our calcu- 
lations is taken from experimental data, and ranged from 
7.1330 A at 100 K to 7.1552 A at 1200 K.^^ For a given 



temperature, a typical run consisted of 10^ MD steps for 
system equilibration, followed by 5 x 10^ steps for the cal- 
culation of ensemble average properties. To have a nearly 
constant precision in the path integral results at differ- 
ent temperatures, we have taken a number of beads, L 
(Trotter number), that scales with the inverse tempera- 
ture such that LT ~ 6000 K. For comparison with the re- 
sults of our PI MD simulations, we have carried out some 
classical MD simulations with the same interatomic in- 
teraction (setting L = \). The quantum simulations were 
performed using a staging transformation for the bead co- 
ordinates. Chains of four Nose-Hoover thermostats were 
coupled to each degree of freedom to generate the canon- 
ical ensemblei2£ To integrate the equations of motion we 
have used the reversible reference system propagator al- 
gorithm (RESPA), which allows one to define different 
time steps for the integration of the fast and slow degrees 
of freedomiiSS, For the evolution of the fast dynamical vari- 
ables, that include the thermostats and harmonic bead 
interactions, we used a time step 6t — At/4, where At 
is the time step associated to the calculation of DF-TB 
forces. The convergence of the total energy as a function 
of At is shown in Fig. ^ A value of At = 0.5 fs is found 
to provide adequate convergence. We have also explored 
the convergence of the simulation as a function of the 
thermostat mass, Q, 



ff3n^ 



(1) 



where /3 = {kBT)~^ is the inverse temperature, and / 
is a scaling factor. The standard deviation of the total 
energy, as derived from a block analysis,^ ~ is displayed in 
Fig. 12 at two different temperatures. We observe that at 
300 K the standard deviation can be reduced by about 
20 % by changing the / factor from a value of 1 to a 
value of 0.2. Taking into account that the standard de- 
viation varies with the number of simulation steps, A/5, 

— 1/2 

as Mg , then, for a given threshold accuracy, a sim- 
ulation run using / = 0.2 at 300 K requires 35 % less 
simulation steps than a run using / = 1. In the simu- 
lations presented below, we have varied the parameter / 
linearly for temperatures between 300 and 1000 K. The 
/ values changed from / = 0.2 (for T < 300 K) to / = 1 
(for T > 1000 K). 



B. Calculation of anharmonic vibrational 
frequencies 

To calculate vibrational frequencies we will employ a 
method based on the linear response (LR) of the system 
to vanishingly small forces applied on the atomic nuclei. 
With this purpose, we consider a LR function, the static 
isothermal susceptibility , that is readily derived from 
PI MD simulations of the equilibrium solid, without hav- 
ing to explicitly impose any external forces during the 
simulation. This approach represents a significant im- 
provement as compared to the standard harmonic (HA) 
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FIG. 2: Standard deviation of the total energy of the diamond 
supercell as a function of the scale factor /, used in Eq. ^ 
to define the thermostat mass, Q. The results correspond to 
simulation lengths of 10^ MD steps at 300 K and 2 x 10^ MD 
steps at 1000 K. The lines are guides to the eye. 



approximation. A sketch of the method is given in the 
following. 

Let us call {Rp} = {2;^^} the set of 3NL Cartesian 
coordinates of the beads forming the ring polymers in the 
simulation cell (i = 1, . . . , 3N;p = 1, . . . , L). We consider 
the set {Xi} of centroid coordinates, with Xi defined as 
the mean value of coordinate i over the corresponding 
polymer: 



X, 



L ^ 



(2) 



Then, the linear response of the quantum system to small 
external forces on the atomic nuclei is given by the sus- 
ceptibility tensor x"^, which can be defined in terms of 
centroid coordinates as^S 



T 



(3) 



is the mass of the nucleus asso- 
/iy = {X,Xj) - {X,){Xj) is the 



where j3 — {kBT)^^, rn 
ciated to coordinate i, 
covariance of the centroid coordinates Xi and Xj, and 
(...) indicates an ensemble average along an MD run. 

The tensor allows us to derive a LR approximation 
to the low-lying excitation energies of the vibrational sys- 
tem, that is applicable even to highly anharmonic situ- 
ations. The LR approximation for the vibrational fre- 
quencies reads 



(4) 



where A„ [n = 1,...,3A'^) are eigenvalues of x^, and 
the LR approximation to the low-lying excitation energy 
of vibrational mode n is given by huJn,LR- More details 
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on the method and illustrations of its ability for predict- 
ing vibrational frequencies of solids and molecules can be 
found elsewhere2fli2ii2Si2i. In connection with the vibra- 
tional modes that actually appear in our calculations, we 
note that the application of periodic boundary conditions 
is physically equivalent to the consideration of lattice vi- 
brations only at the center (k — 0) of the BZ of the 
employed simulation cell. Modes with k ^ violate the 
periodic boundary conditions, because all atomic images 
of an atom have distinct displacements, whose amplitude 
is modulated by both the propagation vector k of the 
vibrational mode and the translational vector of the im- 
age. However, periodic boundary conditions implies that 
all atomic images must display a displacement identical 
to that one of the atom located in the simulation cell, a 
condition that is only met if the propagation vector is k 
= 0,^ 



C. Calculation of one-electron energies 

For the sake of clarity, we use the Schrodinger formu- 
lation to derive the expectation value of electronic ob- 
servables. However, the final result is obtained in a form 
appropriate to the path integral formulation. Within the 
adiabatic Born-Oppenheimer approximation,'^^ the total 
wave function is written as 



*,(r,R) =x»(R)¥'o(r,R) 



(5) 



where (r, R) are the electronic and nuclear coordinates, 
Xi labels the nuclear wave function, and (po represents 
the electronic ground state configuration, which depends 
parametrically on R. Let us call the energy of the state 
^I^i. Then, the canonical partition function, Z, is defined 



(6) 



We consider an electronic observable represented by the 
operator £(r, R), that is a function of the electronic coor- 
dinates and depends parametrically on R. Its canonical 
average, (E), is defined as 



{E) = Z 



dR J dr^*{r, R)f(r, R)*,(r,R) 



. (7) 

Now, we write this equation in an alternative way, bet- 
ter adapted to the path integral formulation. First, the 
function E(R) is defined as the expectation value of the 
operator £(r,R) over the electronic wave function. 



E{R) = J dr(p5(r,R)5(r,R)(po(r,R) • 



(8) 



The second definition is the funtion p(R, R), representing 
the diagonal elements of the normalized canonical density 
matrix for the nuclear coordinates, 



p(R,R) = Z- 



(9) 



Considering the factorization of ^'i(r,R) in Eq. (O, and 
using the last two definitions, we can rewrite the average 
(E) in Eq. © as 



(E) 



dR.p(R, R)£;(R) . 



(10) 



The last equation shows that electronic observables, (E), 
are obtained as ensemble averages over the nuclear con- 
figurations, R, accessible in thermal equilibrium. This 
equation can be readily used in combination with the 
path integral sampling of the density matrix /9(R, R). We 
will apply it to calculate the canonical average, (En), of 
the n-th energy eigenvalue of the electronic Hamiltonian. 
For convenience, let us define the following probability 
density for the one-electron state En, 

Pn{E)^ j dRp{R,R)5{E.n{R)-E), (11) 

where 5 is the Dirac delta function. Pn{E) can be easily 
accumulated during a PI MD simulation. The expecta- 
tion value, {En) , as given by Eq. 1)10(1 . can now be written 
as the first moment of the distribution pn 



{En 



dEpn{E)E. 



(12) 



The direct electronic gap of diamond is derived as 

Eg = {E,) - {E,} , (13) 

where Ec and E^ are the one-electron states associated 
with the bottom of the conduction band and the top of 
the valence band at the reciprocal point k = 0. 



III. RESULTS AND DISCUSSION 

In this section we present the main results derived from 
our PI MD simulations of diamond as a function of tem- 
perature. Whenever possible we will compare the sim- 
ulation results to available experimental data. Results 
concerning the electronic and vibrational properties are 
presented in the next subsections. 



A. Electronic properties 

1. One-electron states 

The top of the valence band, E^, for a diamond crystal 
with the atoms fixed in their crystallographic positions, 
Rmin, is three- fold degenerate. Each of these three one- 
electron levels leads to identical probability densities, pv, 
as defined in Eq. ((TT|l . The function Pv{E), derived from 
our PI MD simulations, is represented in Fig. Oat two 
temperatures, 100 and 1000 K. The probability density 
shows three distinct maxima reflecting that the underly- 
ing electronic state is three-fold degenerate. In each case. 




energy (eV) 

FIG. 3: Probability density functions for the valence band 
top of diamond obtained by PI MD simulations at 100 and 
1000 K. The expectation value of the electronic level, (Ev), 
is shown as a vertical continuous line. The broken line shows 
the valence band top, Ev, for the nuclear configuration Rmin, 
with the atoms fixed at their crystallographic positions and 
the cell parameter fixed at the equilibrium value at 100 K. 
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FIG. 4: Probability density function for the conduction band 
bottom at the reciprocal point k = 0, obtained by PI MD 
simulations at 100 and 1000 K. The expectation value for 
this electronic level, (-Ec), is shown by a continuous vertical 
line. The broken line shows the position of Ec for the static 
nuclear configuration Rmin. 



the continuous vertical line displays the expectation value 
of the valence band top, (-Ei,) , defined as the first moment 
of the distribution p„ [see Eq. (|12f) ]. The temperature ef- 
fect in (Ey) is a shift of 0.03 eV toward higher energies, 
when the temperature increases from 100 K to 1000 K. 
This shift is a consequence of the different magnitude 
of the electron-phonon interaction at both temperatures 
(longer atomic displacements at higher T). The dotted 
line in the figure represents the energy of the top of the 
valence band, E^, for a crystal with atoms fixed at their 
equilibrium positions and cell parameter set to the equi- 
librium value at T = 100 K (a = 3.5665 A). The PI MD 
simulation predicts that at 100 K the top of the valence 
band is shifted by 0.11 eV with respect to the result ob- 
tained for the Rmin configuration. 

In Fig. 21 we present the results corresponding to the 
one-electron state Sc, i.e., the conduction band bottom 
at the reciprocal point k = 0. The energy shifts found for 
as a consequence of the electron-phonon interaction 
are larger and of opposite sign as those encountered for 
the valence band. At 100 K we observe a downwards 
shift of {E^ by about -0.59 eV with respect to the Rmin 
configuration. The temperature shift in (iSc) amounts to 
-0.37 eV between 100 and 1000 K. 



2. Direct electronic gap 

The first direct gap. Eg — {Ec) — (E^), of diamond 
has been portrayed in Fig. [51 as a function of tempera- 
ture. The open circles have been derived from PI MD 



simulations where the thermal expansion of the crystal 
lattice has been taken into account by varying the value 
of the cell parameter a. Open squares are PI MD results 
derived at different temperatures with a cell parameter 
fixed at the equilibrium value at 100 K. We note that 
the effect of the thermal expansion on the electronic gap 
is a slight reduction of the gap. The temperature de- 
pendence of the electronic gap predicted by our PI MD 
simulations is in reasonable agreement with the experi- 
mental results reported in Ref. js^ for diamond Ila up to 
700 K, based on measurements of the complex dielectric 
function by spectroscopic ellipsometry. The slope of the 
simulation results at temperatures above 500 K is larger 
than that of the experimental data. A possible reason for 
this behavior is an overestimation of anharmonic effects 
by the DF-TB potential model. Although we expected 
a reasonable agreement between experimental and the- 
oretical results at low T, the coincidence shown in Fig. 
O for the absolute value of the first direct gap is fortu- 
itous. In fact, the displayed experimental values, derived 
from first-derivative line-shape analysis of the complex 
dielectric function, are shifted by about 0.07 eV toward 
higher energies, in case that they are obtained by second- 
derivative line-shape analysis.-^ 

To quantify the influence of nuclear quantum effects 
on the value of the direct electronic gap of diamond we 
have performed a series of classical MD simulations as 
a function of temperature. The shifts of the energy gap 
obtained in the classical simulations are compared to the 
PI MD results in Fig. The most prominent quantum 
effect appears in the low temperature limit as a conse- 
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FIG. 5: Temperature dependence of the direct electronic gap 
of diamond obtained by our PI MD simulations. Open cir- 
cles were derived from simulations that take into account the 
thermal expansion of the lattice. Open squares correspond to 
simulations where the cell parameter was fixed to the equi- 
librium value at 100 K. The continuous line is the fit to the 
experimental data given in Ref. |^ for diamond Ila. 



FIG. 6: Relative shifts of the direct electronic gap of dia- 
mond. Open circles are results from classical MD simulations, 
while open squares correspond to quantum PI MD simula- 
tions. The solid circle at T = and the continuous line are 
results presented in Ref. based on perturbation theory. 
For both PI MD simulations and perturbation theory, shifts 
are given with respect to the corresponding quantum limit at 
T — 0. The dotted line is a guide to the eye. 



quence of the zero-point vibration. The renormalization 
of Eg amounts to 0.7 eV at T = 0. This value agrees 
well with the perturbational treatment of the electron- 
phonon coupling in Ref. whose result is represented 
by a closed circle in Fig. El We stress that the zero-point 
renormalization of Eg amounts to about 10% of its value. 
The PI MD results show satisfactory agreement with the 
perturbation theory data available up to 700 K,^ The 
main discrepancy found between both sets of results is 
that the slope at temperatures above 500 K is larger for 
PI MD than for perturbation theory. The overestimation 
of anharmonic effects by the DF-TB model is a probable 
explanation for this behavior. Differences between the 
quantum and classical results for Eg are significant in 
the whole studied temperature range. In particular, at 
room temperature the classical result deviates from the 
PI MD data by about 0.45 eV. 



We have calculated the shift of the direct electronic 
gap of diamond as a function of the isotopic mass at 300 
K. The calculated energy gap for ^^C amounts to 7.054 
eV. For ^^C, this gap increases to 7.081 eV at the same 
temperature. The isotopic effect of 27 meV at 300 K 
is in reasonable agreement to the value reported in Ref. 
of 22 meV, based on perturbation theory in the zero 
temperature limit. 



B. Vibrational Properties 

1. Vibrational energy 

The vibrational energy of the simulation cell of dia- 
mond, as derived by the PI MD simulations, is presented 
as a function of temperature in Fig. [7| (solid circles). The 
zero of the energy scale corresponds to a diamond crys- 
tal with fixed atoms and with the cell parameter fixed 
at the equilibrium value at 100 K. The thermal occupa- 
tion of excited vibrational states is evident in Fig. Qby 
the increase of the vibrational energy with temperature. 
To quantify the anharmonic effect on the vibrational en- 
ergy of the crystal, we have plotted in Fig. [7|the har- 
monic vibrational energy (solid diamonds). The set of 
harmonic frequencies has been derived by diagonalizing 
the dynamic matrix in a simulation cell with the exper- 
imental equilibrium lattice parameter at each tempera- 
ture. At the lowest studied temperature the harmonic 
result deviates from the PI MD value by about 0.1 eV, 
which amounts to about 0.8 % of the total vibrational 
energy. This error of the harmonic approximation is a 
consequence of the anharmonicity of the phonon vibra- 
tions. 

The set of anharmonic vibrational frequencies, uJn.LR, 
derived by our linear response approach, are expected 
to represent an improved description of the vibrational 
problem. Thus, we have recalculated the harmonic vi- 
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FIG. 7: The vibrational energy of the simulation cell of dia- 
mond as a function of temperature is shown by open circles. 
The zero of energy corresponds to the crystal with fixed atoms 
and cell parameter a = 3.5665 A. The harmonic vibrational 
energy derived from the set of harmonic and linear response 
frequencies is displayed as open diamonds and squares, re- 
spectively. The lines are guides to the eye. 



brational energy by considering LOn,LR as a set of renor- 
malized phonon frequencies. The result is shown as open 
squares in Fig. [7| Most of the error of the harmonic 
approximation is corrected by the LR frequencies. At 
low temperatures the absolute error of the improved es- 
timation of the vibrational energy amounts to 0.03 eV, 
i.e., about 0.2 % of the PI MD result. Our conclusion 
from this comparison is that the set of LR frequencies 
provides a consistent description of anharmonic effects 
in the employed DF-TB model, in Hnc with previous 
results of vibrational properties on molecular and solid 
state systems ^^Si^SiS In the following, we focus on the 
study of two particular LR phonons that can be com- 
pared to available experimental data. 



2. Optical phonon at k = 

The highest energy phonon of the diamond crystal is 
the optical phonon at the center of the BZ (k = 0). At 
200 K, the LR wavenumber of this phonon amounts to 
1396 ± 5 cm^^. The harmonic result, obtained for the 
equilibrium cell parameter at this temperature, is 1407 
cm^^. The difference between the LR and harmonic re- 
sult is a consequence of the anharmonicity of the inter- 
atomic potential that induces a slight softening in the 
phonon frequency. The optical phonon wavenumber of 
diamond in the zero temperature limit is found in the 
first-order Raman spectrum at 1335 cm~^ii The DF-TB 
potential overestimates the experimental wavenumber by 
about 61 cm~^. This error is lower than that found in 
other tight-binding parametrizationsi^SiS 
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FIG. 8: Temperature shift of the optical phonon at F (k = 0) 
of diamond. The results of the LR and harmonic approxima- 
tions are compared with a fit to the experimental data of Ref. 
0- The extrapolated LR value at T = has been employed as 
zero of the vertical axis for both the LR and HA results. The 
zero of the experimental data corresponds to the extrapolated 
experimental wavenumber at T = 0. 



The relative shift of the optical phonon in diamond is 
presented as a function of temperature in Fig. |S1 The 
comparison between the harmonic and LR results shows 
that anharmonic effects lead to a softening of the phonon 
mode in the studied temperature range up to 1200 K. 
The LR results deviate from the experimental data at 
temperatures above 500 K. This deviation shows that the 
anharmonicity of the optical phonon is overestimated by 
the DF-TB Hamiltonian at those temperatures. We have 
previously commented on the enhanced anharmonicity of 
the DF-TB model in relation to the decrease of the direct 
electronic gap at temperatures above 500 K. 



3. Elastic constant C44 

The lowest energy phonon, wi, determined by either 
the LR or HA approximations is 12-fold degenerated in 
the employed simulation cell. The wave vector, ki, of 
this phonon state has been identified as the midpoint 
between F and X points along the A [100] direction in re- 
ciprocal space, with coordinates (7r/a,0, 0). This means 
that the phonon oji corresponds to the transverse acous- 
tic (TA) branch. The identification of the wave vector 
of oji has been possible by solving the dynamical ma- 
trix at some selected k points of the BZ of the primitive 
unit cell of diamond. The phonon velocity, vtAi along 
the TA branch, l^taC^ta), is defined as the slope of the 
dispersion branch at the origin. 



ujta ^1 
VTA = hm « c— 



(14) 



In this equation, ki = ir/a is the modulus of the wave 
vector ki associated to the phonon uji. The constant 
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FIG. 9: Temperature shift of the elastic constant C44 of di- 
amond. The results of the LR and harmonic approximations 
are compared to a fit to the experimental data of Ref. IsM 
The extrapolated LR value at T = has been employed as 
zero of the vertical axis for both the LR and HA results. The 
zero of the experimental data corresponds to the extrapolated 
experimental constant at T = 0. 

c — 1.037 corrects, for the case of the HA approximation, 
the finite difference error encountered by using the finite 
vector ki to calculate the slope at the origin. The elastic 
constant C44 can be derived from vta by the relation 

C44 = PVXA : (15) 

where p is the density of the diamond crystal. 

Using the harmonic value of oji^ha at 200 K we obtain 
with Eqs. idJ) and a value for C44 of 551 GPa. The 
estimation of the elastic constant using the LR wavenum- 
ber, LOi^LR, in Eq- (|14|l is of 545 GPa. The experimental 
result derived from Brillouin scattering amounts to 576 
GPai^ The shift of the elastic constant C44 with temper- 
ature is plotted in Fig. O The comparison of the LR and 
harmonic results shows that anharmonic effects cause a 
reduction of the value of the elastic constant C44 in the 
whole studied temperature range. The comparison to the 
experimental data shows that the DF-TB model gives a 
reasonable prediction of the temperature shift of the C44 
elastic constant, even though there appears a systematic 
underestimation of this shift. 



have chosen a simplified electronic Hamiltonian to de- 
velop the algorithms required for the simulation of solid 
state systems, but this limitation should be eliminated in 
the future by the implementation of improved electronic 
structure methods. 

The temperature and isotopic dependence of the first 
direct gap of diamond predicted by our PI MD simulation 
shows good agreement with the available experimental 
results, based on spectroscopic ellipsometry^ and theo- 
retical results, based on perturbation theoryi^ Thus, the 
employed simulation model has demonstrated its capa- 
bility to realistically describe electronic properties that 
are determined by electron-phonon interactions. The ef- 
fect of the zero point vibrations of the lattice phonons of 
diamond in its first direct gap is a reduction of the gap by 
about 0.7 eV. This effect is so large that any theoretical 
approach aiming at a quantitative determination of the 
electronic gap of diamond can not be based only on an 
improved solution of the many-body electronic problem, 
but it should also include the treatment of the electron- 
phonon coupling. 

Anharmonic effects in the lattice vibrations have been 
derived by a linear response approach based on the path 
integral formulation. This approach allows us to de- 
rive anharmonic vibrational frequencies from the study of 
spatial correlations in the displacements of the vibrating 
nuclei in the simulation cell. Anharmonic effects are re- 
sponsible for a reduction of the vibrational energy of the 
solid of about 1 %, with respect to the result predicted by 
a harmonic approximation. The temperature shift of the 
optical phonon at k = is larger that the experimental 
result determined by Raman spectroscopy.^ We consider 
that this limitation is a consequence of the parametriza- 
tion of the employed DF-TB one-electron Hamiltonian^SS 
that overestimates the anharmonicity of the highest fre- 
quency phonon of the lattice. Better agreement is found 
in the comparison of the elastic constant C44 derived from 
the simulations with the experimental values obtained by 
Brillouin scattering 1^ 

We plan to extend our simulations to more complex 
systems like hydrogen impurities in diamond, ■^^ where the 
presence of light impurities should strengthen further the 
influence of quantum nuclear effects in the electronic and 
vibrational properties of the lattice. 



IV. CONCLUSIONS 

The simulations presented here for diamond are based 
on the treatment of electrons and nuclei as quantum par- 
ticles in the framework of the Born-Oppenheimer approx- 
imation. The use of the path integral formulation for the 
atomic nuclei allows us to obtain the vibrational and elec- 
tronic properties of the solid at finite temperatures. We 
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